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DEVELOP MINIMUM THRUSTOR CONTROL LAWS AND SELECT ORBITS 
FOR A GEODESY DRAG -FREE SATELLITE 


INTRODUCTION 

The research proposed in Ref, 1 and subsequently made a part of the Work 
Statement included two specific areas: Research in control laws for pulse 

plasma systems, and mission analysis. These are discussed separately below. 
We have summarized principal results but reference the progress letters [2] 
for many details. In the discussion of the pulse plasma control laws, 
we have also included some work that has subsequently been done under Ref. 

3. 


A. PULSE PLASMA CONTROL LAWS 

The original motivation for studying control laws for pulse plasma 
systems was based on the improved life characteristics possible with pulse 
plasma jets. These pulse plasma units are relatively massive compared with 
cold gas thrustors. As a result, therefore, significant mass savings can 
be achieved by minimizing the number of thrustors. The control laws, there- 
fore, were developed for thrust available from two thrustors only. In a 
spinning satellite, these thrustors are sufficient to completely control 
the vehicle as long as the spin rate is sufficiently high for a given level 
of external disturbance. The ^ rustnrs are canted so that a component 
each is along the + and - spin axis. The other component of each thruster 
acts in the radial direction. It is sufficient to analyze the behavior in 
the plane of spin assuming a single thrustor. The motion along the spin 
axis is equivalent to a drag-free control system of the type that was flown 
on the successful DISCOS in 1972 September, 

It became clear that this type of control system would be adaptive to 
a cold gas system also as it can result in considerable simplification. 
Perhaps more importantly, by having more than one thrustor in the rotating 
plane, a completely redundant system can be designed with no penalty in 



hardware over what is normally used lor three-axis control. We have there- 
fore' continued our studies beyond the scope of the present contract with 
innovation for the cold gas thrust system to improve the reliability of the 
system design being studied under Ref. 3. The early results that were 
reported in Progress Letters 4 and 5 [of Ref. 2] form the basis on which 
the more recent research has been performed. 

In summary, we have been working with a control law based on neutral- 
izing the velocity with respect to an inertial reference and providing a 
small component of velocity which drives the proof mass toward the center 
of the cavity. This algorithm is based on the fact that action can only 
be taken when the single thrustor is correctly aligned and typically pro- 
vides impulses with approximately a one-spin- period frequency. 

Figure 1 illustrates the dynamical model of the system in which we have 

the satellite stably spinning in a locally inertial environment. The motion 

in the spin (z) direction is assumed uncoupled and controlled. The spin 

plane body axes x^ y fe are centered at the satellite mass center (c.m.). 

Initially the position sensor null is also located here. The control action 

is applied by a component of force from a jet fixed to the satellite along the 

-x direction and directed toward the origin. The position of the proof mass 
b _ 

is determined by its measured and y^ offsets. In addition, an auxiliary 

frame N, (x , y ) , is fixed at the satellite mass center, initially co- 
incident with x b , at some time t Q and maintaining the same orientation 

relative to some inert ially fixed reference. Because of the symmetry of the 
dynamical model, we treat the motions of the system as if the satellite mass 
center were inertially fixed and the control were applied to the proof mass 

in the -£ direction. The magnitude of the force is of course scaled by 
b 

the mass involved. Finally, we assume that the velocity of the proof mass 
relative to the satellite c.m. in the N frame and resolved along the x^ t $ b 
axes is available from a suitable estimator. 

The control force is analyzed using two models. First, the control force 
(AV ) may be specified proportionally above a particular minimum value (AV ). 
Also, a pulsed control force is considered in which an integral number of pulses 
of size AV ni may be specified. The applicability of these models is deter- 
mined by the size of AV™ and the disturbance environment. 
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Initially, the analysis considers two types of disturbances • First, 
an initial proof mass position and velocity is considered. Also, a constant 
external force (e,g., drag) is assumed in the -x^ direction producing a 
"force” on the proof mass in the + direction. 

The control logic considered asks that the velocity of the proof mass 
after control action in the N frame be directed toward the satellite c.m, 
and be proportional to the distance from the c.m. 

Define V + = velocity of proof mass in N after control 
r “ (x^, ^b^ = P osit ^ on vector of proof mass 

to = satellite angular rate (constant) 
k = proportionality constant. 

Then, the control law would require 



k 

— 7— ** 
JC/U - 


( 1 ) 


where k measures the fraction of the offset distance travelled by the proof 
mass in a half satellite- revolution. If we define 


V~ = velocity of proof mass in N just prior to control 
AV^ = required control action 


then 



Again, defining 

a/ £ effect of applied control action on velocity of proof mass, 
if we assume the proportional- above -a-rainimum control, then from (1): 
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If we assume the integer-pulse control, several logic schemes are avail 
able. The first one considered specifies the control applied to be the great- 
est integer dividend of | AV^I by AV m : 


Cl II: 



AV R « x > 0 

— b ~ 

AV R - x, < 0 

— b 


where G(x/y) = greatest integer dividend of x by y. If we define 



and assume that f is randomly distributed between 0 and 1, then 

E( f ) = I . 


This average control deficiency manifests itself as a c.m. offset and to com- 
pensate for this effect, a third control logic might be specified: 


CL III: 




AV R * x > 0 
— b — 


AV 


.m 



A. 



< 0 


In this case E(f) = 0 and the average control force is correct. 


It was assumed in the analysis that the time history of any control 
application was much shorter than the spin rate. Thus, the control action 
consisted of an impulsive AV applied of magnitude [aV^. Since the 
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orientation of the jet is ‘fixed to the body, the direction of AV^ is always 
directed in the -x direction. Triggering information is derived from the 
calculation of the component of AV R along the y b axis. When this com- 
ponent is zero, the jet is aligned for firing. The general vector component 
control logic is as follows: 

f v * b + Tfo x b Av ^ b > 0 

a/ = I 

X b [ 0 AV^ < 0 

AV^ s 0 



’CL 




where g is one of the control logic functions specified above. 

C Li 

A hybrid simulation of the plant and controller was assembled on our 

HP 2114/EAI-TR48 system to determine the performance of the control laws 

above. Figure 2 depicts this simulation. The dynamics of the proof mass 

/ 2 

motion were simulated in inertial coordinates, i.e,, two uncoupled 1/s 
plants. The output was represented by an inertial position and velocity vector 


« 

x 

n 



A rotation vector was generated by an oscillator to 
ing signal for a coordinate transformation from inertial 
coordinates, e.g,, 


produce the driv- 
to body-fixed 
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FIG. 2 SINGLE THRUSTOR HYBRID SIMULATION 
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The position and velocity vectors in body coordinates were then 
input to a sensor model and an analog controller to form the error vector, 


- 


k n 

e 

x b 


-v ' ^ x b 

b 



- k 

e 



X. 


y jt/w b 

l *L 


L b J 


which corresponds to the required velocity. 
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Thus far the sensor model has been linear with only a null bias offset 
from the body axis origin (c.m. location). 

The e signal is compared to a threshold (nominally zero) and a 
yb 

threshold crossing produces a triggering signal to the digital computer. 

The primary function of the digital computer is to simulate switching 
control logic. Numerical calculations for the control function are minimal. 
This is done in consideration for future analog or hard-wired digital cir- 
cuitry which would be resident in an actual vehicle controller. The 




signal is sensed by the computer and using the control logic (e.g,, CLI), 
the appropriate number of thrustor pulses is calculated. Then a pulse train 
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is applied to the inertial plant according to a force vector decomposition 


into inertial coordinates driven by the rotation oscillator 


- - 


“ 
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F C 

X 
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C 

= |f c | . 

cos 
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sin 

tot 

_ n _ 



- 


where 



> 


the control force vector. 


A pulsed plasma thrustor with a microsec pulsewidth, 40 plb-sec impulse 
size, and 2,5 pulses/sec repetition rate was used as a nominal actuator 
plant. 


Aside from its primary control function, the digital computer was 
used as a monitoring system for various performance calculations such as fuel 
consumption and stability. The ability to perform experiments, take data, 
and calculate results automatically, provide the motivation for the hybrid 
simulation adopted. 

Several qualitative results of the simulation were almost immediately 
obvious. These are listed below: 

1, A finite repetition rate tends to spread the direction of applied 
force. This effect tends to destabilize the system for large 
initial conditions and values of k. 

2, An adaptive threshold or threshold with memory can stabilize 
the system for large initial conditions and small repetition 
rates without a large external force. An external force seems 
to quickly drive such a system unstable. 

3, As k -> 1.00, the system destabilizes in the presence of delays. 
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4. A dead zone on the position (x^ y^) vector tends to drasti- 
cally reduce the limit cycle size and fuel consumption. It 
adds stability with long pulse separations. The center of 
mass (c.m.) hangoff with an external force is not significantly 
affected. 

5. A sensor null bias has a large destabilizing effect. This phe- 
n omen a is not well understood yet and a stability analysis 

has been initiated, 

6. There is a hangoff of the average proof mass position from the 
sensor null in the presence of an external force as predicted 
by the conservation of momentum. The trajectories tend toward 
one-dimensional as can be seen in Fig. 3. For a smaller 
external force, the effect is not as pronounced. See Fig. 4. 


The fuel consumption of the system was calculated at various values of 
external force and initial conditions. The results compared quite favorably 
with the digital computer simulation and are shown in Fig. 5. At higher values 
of the external force, the number of thrustor firings becomes high and the 
spreading effect of a finite repetition becomes pronounced. This causes an 
impulse in the direction transverse to the drag force. This impulse must 
be taken out if the cross-wise momentum is to be conserved. Using this fact, 
a correction (to the fuel consumption) for a finite repetition rate can be 
derived. (see sketch below.) 
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SFC (PULSES A REV) 



FIG, 5 FUEL CONSUMPTION FOR VARIOUS PULSE WIDTHS SIMULATION RESULTS 


- 12 - 



Define: Y 


to 


0 

F C 



AV 


m 


A 

A 

A 

A 

A 

A 

A 

A 

A 


repetition rate (pulscs/sec) 
satellite spin rate 
Y/w s rate factor 
desired control force 
resultant control force 
transverse velocity component 
external acceleration 
number of firings 

minimum velocity change (l impulse) . 


Then, 


a = n/2£ 

AV T « n . AV m . a 
from intrack conservation of momentum: 


or 



But, 


where 


m ~ added fuel consumption (pulses/J rev) 

2AV 

= l/p(AV e ) 2 
AV e £ aVav” . 
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Finally, 


SFC . = AV e + — AV e (pulses/^ rev) . 

mm p 

These values are shown in Fig, 5 and seem to correlate very well with the 
simulation results. 

We conclude that single thrustor control in the spin plane is (1) 
stable for the conditions studied, (2) yields fuel consumption approxi- 
mately equal to the drag except for very low drag levels or spin rates, 
and (3) is practical to mechanize. 
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B. MISSION ANALYSIS 


INTRODUCTION 


Radio altimetry measurement of the ocean surface from a rather low 
altitude satellite with an accurately known orbit makes it almost mandatory 
that this orbit be both drag-free and approximately circular. The accurate 
determination of this low altitude orbit (i.e., determination of the geo- 
centric radius to better than 1 meter) which requires, of course, calibra- 
tion of tesseral harmonic perturbations which are not accurately known 
a priori, can be accomplished either with range measurements between a 
pair of drag-free satellites in the same orbit, or with range measurements 
from a much higher satellite, possibly also drag-free, on which tesseral 
perturbations are much smaller and, hopefully, removable. 


ANALYTIC .DESCRIPTION OF THE LOW ORBIT 

This section deals with an analytic description of a near circular, 
ideally drag-free, orbit. This description is especially suited to the 
satellite pair ranging method of orbit determination, but it is useful 
also in the context of orbit determination by range measurements from a 
high altitude satellite, especially if the high satellite remains in essentially 
the same plane as the low satellite, e.g., if both orbits are perpendicular 
to both the equator and the ecliptic planes. The orbit parameters employed 
are not the conventional ones, but are related to parameters used in 
earlier research at Stanford University [4, 5], The effect of tesseral 
harmonics is described directly in terms of position perturbations 
rather than orbit element perturbations, and these are expressed in terms 
of certain coefficients which are combinations of the conventional C £,m’ 

S j£m' 

It is assumed that all tesseral resonances (up to fairly high order) 
are avoided. This assumption about the choice of the low altitude orbit 
guarantees that no small unknown terms in the expansion of the geopotential 
can have an unduly large effect on the orbit. An analytic description, 
which can be good, after fitting, to better than 1 meter is then quite 
manageable and is represented by the following equations. 
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— = -fc cos 0 - ri sin 
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1 / R ®\ . 2 . 

9 + 4 J 2 \t' 51 " 1 


in i cos 20 + (9 r 


R ' 


1 / ' 0 | ^ 

50 = 2£ sin 0 + 2r) cos 0 + - J l — / sin i sin 20 4- ( 9, 

O Z 3. 


where a is the mean radius ( not the average osculating semi-major 
axis), and ( §, rj) is a "mean eccentricity" vector. Here, 0 is a mean 
orbit position given by 


0 = nt + constant , 

where n (the "nodal mean motion") is given in terms of a by 


(!js) f 1 + J (V) 

(| 

- H + *) 

V 3 / l 2' a ' 

V 4 

8 / 2 ) 


It will be useful to assume that i is the inclination of a "best fitting" 
mean plane regressing at a constant rate ( - &) so that the earth's oblate- 
ness gives out-of-plane excursions only of order J^| or . This i 

is not the average osculating inclination. With this definition of the 
mean plane, out-of-plane excursions are negligible in their contributions 

to the range between satellites in the same or nearly the same plane. 

2 

The terms 0 ^ in the formula for n include terms such as J„, J., J , 

2 2 ; 4 J 6' 

as well as terms from luni-solar perturbations. All these terms 

can be obtained from Dasenbrock [4], after finding the relation between 

-2 , 

the present mean radius a and his parameter h f\ i. The terms 0^ 
in 8r, 50 include terms of second order jointly in J , £, and Tj, as 
well as terms from J^( i «= 2), C. , (m ^ l), and from luni-solar perturba- 

tions. The tesseral harmonic perturbations may be expressed in the form 
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where X is the longitude of the ascending node, relative to Greenwich, 
N 


Here, 


km 


cosj 

u 

and 

B, \ 

<• , 

cos 

- 
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[ km 


km 

sin ^ 
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km 


are known linear combinations of the usual C , S^, with & - 

jk|, j kj + 2, | k | +4 ••• [see Morrison, 5], 

The mean eccentricity vector (|, Tj) undergoes long period changes: 


60 


in 

dO 


* 

**u) rj + k 


-* 


= to 


where 


to = J 


£)< 


2 - 


5 2 2 

— sin i) + C?(j , J. , •»«) + luni- solar terms 

2 2 4 b, 


and 


k = linear combination of J , J j ... 

o' 5 ; 7 ' * 

all of which are obtainable from Dasenbrock [4], 
THE SATELLITE PAIR RANGING METHOD 


Figure G illustrates the relative geometry of the satellite pair in orbit 
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"MEAN" ORBITAL PATH 



FIG 


6 GEOMETRY OF THE SATELLITE RANGING EXPERIMENT 



The fluctuations in the measured range R are expressible as 


SR = ( &r^ 4- Sr^) sil1 \ + a ( ^9 “ cos *2 + ^ 2 * 


It is assumed that both satellites are launched together and then drift apart 
very slowly so that the separation angle OL increases from zero up to a 
maximum permissible separation OL , at which the radio wave passes 
within 70 km of the earth's surface. For greater separations, the atmos- 
pheric refraction will introduce pathlength errors of the order of 1 meter. 

As soon as a reaches OL . one of the satellites is given a small 

max 

AV either to stabilize CL or to reduce it again towards lower values. 

Prior to reaching OL , the continuous measurement of R provides 

max 7 

enough information to determine not only both eccentricity vectors but also 

all of the coefficients A and B up to a high order, with the 

km km 

exception of the coefficients which do not affect R, The AV 

application will introduce a further small change in one of the eccen- 
tricity vectors, which, however, is rapidly estimated. 


REFINEMENT OF THE EARTH'S GRAVITATIONAL CONSTANT 

A side benefit, which is possible from the satellite pair ranging 
experiment, is the refinement of the earth’s gravitational constant (p ), 
although, as will appear, this refinement improves markedly with an in- 
crease in the altitude of the satellite pair. 


Present uncertainty in p is about three parts in 10 , although 
tracking data from Mariner 9 may have reduced this uncertainty by a 
factor of 2. Now the nodal period of the satellite pair is obtainable to 
very high accuracy from ground tracking. This permits the determination 

7 

of p to one part in 10 if the mean radius a of either satellite can 

7 

be determined to better than one part in 3 X 10 . It may be assumed that 

the correction terms in the expression for n would be known to better 

4 

than one part in 10 from current values of •••, together with 

ground tracking determination of i (and a). 
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To determine a with the required accuracy; consider the measured 
range with the short period fluctuations averaged out 

R av. “ 2a Sin f [ X + ® (J 2 > + ^ 

where CL } the mean separation angle, is expressible as 

a = at . 

The accuracy with which a itself is estimated depends strongly on the 
maximum useful CL, since for small CC, R is essentially the product 
( aGl) t which is insufficient to estimate a and CC separately. This 
accuracy is indicated in Figure 7 , showing the bias in estimating a 
due to a constant bias AR in the range measurement, plotted vs useful 

a. 

It has so far been assumed that both satellites are "drag-free", 
so that the intrack position contains no acceleration due to atmospheric 
drag. Recent orbital success with a drag-free satellite indicates that 
the mass attraction of the satellite on the proof mass produces an in- 
track acceleration of order of magnitude 10 ^g. If we assume gravity- 
gradient stabilized satellites, these accelerations may be assumed essen- 
tially constant, but unknown a priori. The estimation of the relative 
acceleration, together with that of a and a, slightly degrades the 
accuracy of the estimate of a (or more correctly, of the mean a during 
the interval of observation). This is included in Figure 7, However, 
small fluctuations in the intrack accelerations about their mean values 
will further degrade the accuracy of estimation of a. 

The mean intrack accelerations can be further reduced if the satel- 
lites are spun about axes perpendicular to their orbit plane. If the 
spin period is comparable with the limit-cycle period of the drag-free 
translation control system, the mass attraction may be expected to fluc- 
tuate, but the intrack relative acceleration A should now have essen- 
tially zero mean and a fluctuation which may be crudely modeled as a 
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ORBITAL ALTITUDE ABOVE EARTH'S SURFACE (km) 


FIG. 7 BIAS IN THE ESTIMATE OF SEMI -MAJOR AXIS DUE TO CONSTANT BIASES 
IN RANGE MEASUREMENT 


~ 21 - 


random process with standard deviation of order of magnitude 10 g 

and a correlation time t. of order of the spin period. 

A 

The effect of such a random relative acceleration on the separation 
Of and hence on the estimate of a is summarized in Figure 8, 

From Figure 7 we see that it will suffice to reduce the bias AR 
in range measurement to 10 cm, if the orbital altitude exceeds 1000 km. 

From Figure 8, if we assume, for example, an altitude of 1000 km 

and a separation rate of 6 x 10 rad/sec = 3°/day, then using 

a ~ 10 11 g, t. = 1 sec, cu/a = 0,3(10 7 ). 

A Ad 

It should be added that Figure 8 was calibrated for the average fit 
of a noise-free model to data obtainable with a noisy intrack accelera- 
tion, It is for this reason that a slower separation rate © would 
actually increase the error of the estimation. An optimal filtering 
scheme would not have had this defect. 
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FIG. 8 RMS ERROR IN THE ESTIMATE OF THE SEMI -MAJOR AXIS DUE 
TO SMALL RANDOM INTRACK ACCELERATION. 
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